A PRECISION RECURSIVE ESTIMATE FOR 
EPHEMERIS REFINEMENT (PREFER) 


Bruce Gibbs 

Business and Technological Systems, Inc. 


ABSTRACT 

PREFER is a fil ter/ smoother program for orbit determination which is used to refine the ephemer- 
ides produced by a batch least squares program (e.g., CELEST). PREFER requires, as input, a file 
containing the nominal satellite ephemerides and the state transition matrices as generated by 
CELEST. PREFER interpolates from this file at the times given on the Measurement Data file and 
processes the measurements in the Kalman Filter to estimate the corrections to the nominal trajec- 
tory. The filter state also includes other parameters which have an effect upon the orbit determina- 
tion (e.g., drag, perturbing gravitational accelerations, thrust, measurement biases and refraction 
parameters, etc.).. Because PREFER is estimating the corrections to the nominal values, all partials 
are evaluated about the nominal trajectory and the filter is linear (not extended). 

The measurement data types which PREFER can process include ground range, range difference 
and Doppler measurements, GPSPAC pseudorange and pseudodelta-range measurements, NAVPAC 
range difference measurements and altimeter measurements. A GPS Trajectory file supplies the 
ephemerides of the GPS satellites which are required to process the GPSPAC or NAVPAC measure- 
ments. A unique feature of the program is the capability to estimate hundreds of pass-disposable, 
measurement biases while using storage and computation for only a few biases. 

After running the Kalman filter forward to the end of the Measurement Data file, PREFER performs 
optimal smoothing. A file created by the Kalman filter is read backward in time and the smoothed 
estimates are obtained by using the recursive formulation of Rauch-Tung-Striebal. 

The combination of a Kalman filter and a smoother should result in greatly improved estimates of 
satellite ephemerides as compared to the batch estimation. Batch estimation is subject to errors 
because of errors in the dynamic models (e.g., gravitational). A filter/smoother which properly 
accounts for dynamic (state) noise should weight the data optimally and reduce the estimation 
errors. Smoothing will produce better estimates (in the middle of the data span) than just a 'forward 
filter because past and future data is used to estimate the state at each point in time (a filter uses 
only past data). Smoothing also tends to average out any dynamic modeling errors which remain. 


PREFER s capability for improving orbit determination has been demonstrated on simulated data 
which contained significant modeling errors. The nominal trajectory had errors as large as 53 
meters and the GPS trajectory file had peak errors of 12 meters. However, the PREFER smoother 
estimate was usually accurate to 3 meters with peak errors of 8 meters. Even during data gaps, the 
smoothed radial error was always less than 6 meters. 
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INTRODUCTION 


A recursive filter/smoother orbit determination program has been 
developed to refine the ephemerides produced by a batch orbit deter- 
mination program (e.g., CELEST, GEODYN). PREFER can handle a variety 
of ground and satellite-to-satell i te tracking .types as well as satel- 
lite altimetry. It has been tested on simulated data which contained 
significant modeling errors and the results clearly demonstrate the 
superiority of the program compared to batch estimation. 


Input 


The input to the program consists of four files and card input. 

A file containing the nominal (batch estimate) host satellite ephemerides 
and the 6 by 6 state transition matrix (from epoch osculating elements 
to current cartesian elements) is interpolated at the times given on 
the measurement data file. A GPS trajectory file supplies the ephemerides 
of the GPS satellites which are required to process the GPSPAC or NAVPAC 
measurements. A sun/moon file supplies the data which is used in the 
earth motion model (for ground based measurements). The card input to 
the program specifies run constants (e.g., time intervals) and a priori 
standard deviations, state noise spectral densities, time constants, etc. 



Measurement Types 


PREFER can process the following types of measurements . 

• Ground Tracking 

Satellite to ground range 

Ground laser range 

Satellite to ground range difference 

Ground Doppler 

• Satellite-to-Satellite 

GPS pseudo range and pseudo delta range 

NAVPAC range difference 

• Altimetry 

Range to center of earth. 

Provisions have been made for handling 50 ground stations and 24 
GPS satellites but only 4 ground stations and 15 GPS satellites can be 
simultaneously observable. This restriction is imposed because of a 
limitation on the total number of states. Since station position 
errors, measurement biases, refraction parameters, GPS position errors 
and timing biases can all be estimated, the state vector could become 
unwieldly. PREFER has the capability to estimate all these parameters 
while using storage and computation for only those parameters which are 
simultaneously observable. This is discussed in later sections. Thus, 
the limitation is on the number of simultaneously observable stations 
and GPS satellites. As a practical matter, this limitation is not very 
restricting since it is unlikely that more than four ground stations 
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would see a low altitude satellite. Furthermore, simulations have 
shown that for the 24 satellite GPS system, no more than 15 GPS satel- 
lites would be observable to a low altitude satellite (without encounter- 
ing severe refraction problems). 

The altimetry measurements are assumed to have been preprocessed 
with a nominal geoid model so that they are treated as a range to the 
center of the earth. 

Dynamics 

A list of the dynamic parameters which PREFER can estimate is given 
below: 

1 Satellite semimajor axis at epoch 

2 Satellite eccentricity x sin (argument of perigee) at epoch 

3 Satellite eccentricity x oos (argument of perigee) at epoch 

4 Satellite inclination at epoch 

5 Satellite mean anomaly plus argument of perigee at epoch 

6 Satellite right ascension of ascending node at epoch 

7 Satellite drag coefficient 

8 Perturbing gravitational acceleration (vertical) 

9 Perturbing gravitational acceleration (cross-track) 

10 Perturbing gravitational acceleration (along-track) 

11 Acceleration of 1st thrust segment (vertical) 

12 Acceleration of 1st thrust segment (cross-track) 

13 Acceleration of 1st thrust segment (along-track) 

14 Acceleration of 2nd thrust segment (vertical) 

15 Acceleration of 2nd thrust segment (cross-track) 

16 Acceleration of 2nd thrust segment (along-track) 

17 Host satellite clock timing error 

18 Host satellite clock drift rate 

19 Altimeter bias. 
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The first 6 are epoch osculating elements. The drag coefficient, 
perturbing gravitational accelerations, host clock drift rate and 
altimetry geoid error (bias) are all assumed to be independent, first 
order Markov processes. This may not be strictly true but it is a 
reasonable approximation. The thrust accelerations are assumed to be 
constant since the thrust durations will be relatively short. 

The state transition matrix for the entire system of dynamic para- 
meters and measurement related biases is: 


4 > = 



0 

I 


where is: 
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The upper left 6x6 partition of is an identity matrix when 

$ is being used to perform the time update on the state vector. However, 
when individual measurements are being processed, the satellite position 
and velocity in cartesian coordinates at the measurement time must be 
known. The nominal position and velocity and the transition matrix from 
epoch osculating to cartesian elements are obtained by interpolation 
from the host trajectory file. The filter state (which includes the 
estimated correction to the epoch osculating elements) is multiplied by 
(ji-j to obtain the estimated correction to the nominal cartesian elements. 

The upper right partition of «j>^ (i.e., the transition from , 

gravitational accelerations and thrust to cartesian elements) is 
obtained as an iterated, second order Taylor series. Since the integra- 
tion time interval will be relatively short (less than 120 seconds) and 
state noise is included in the formulation, a highly accurate integra- 
tion method is not required. 

The state noise covariance matrix (required by the filter) is 
obtained by Taylor series integration of the input spectral density 
matrix. 

Kalman Filter 


Measurements are processed in a Kalman filter to estimate the 
corrections to the nominal trajectory. All partial derivatives are 
evaluated about the nominal trajectory and thus the filter is linear 
(not extended). 

Since the program was intended to process many thousands of 
measurements, the execution time would have been excessive if the 
Kalman equations were evaluated for each measurement. Therefore, the 
measurements are processed in small "mini-batches" (typically 120 
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seconds), during which time, the dynamics are assumed to be deterministic. 
Only when proceeding from the epoch of one mini-batch to the next is 
state noise included in the covariance equations. The term "mini-batch" 
is intended to indicate the lack of state noise rather than the method 
of processing since the estimation algorithm is actually the recursive 
U-D algorithm of Bierman [1], 

A unique feature of PREFER is the capability to estimate hundreds 
of pass-disposable measurement-related biases while using storage and 
computation for only a few. As measurement data from new stations or 
GPS satellites is processed, the state vector and covariance matrix are 
augmented with the a ■priori information for the new measurement para- 
meters. When the station or GPS satellites are no longer visible to the 
host satellite, the parameters are dropped from the state vector and 
covariance matrix. These parameters can be deleted from the filter 
state since they will no longer have an influence on the estimation of 
"common" parameters (dynamic and other measurement related biases). 
However, the deletion of parameters from the filter state does complicate 
smoothing since the lost information must be reconstructed later. This 
is discussed in another section. 

It should be noted that these hundreds of measurement related 
parameters are probably not observable in a statistical sense, i.e., 
a priori information is required to make the covariance matrix full 
rank. These parameters are included in the filter state primarily to 
assure proper weighting of the measurement data. 

Figure 1 is a flow chart of the FILTER subroutine. This routine 
is called once for each mini-batch of data. The flow chart shows the 
sequence of events required to perform the time update, write information 
on the disk for smoothing, process data with the U-D algorithm and 
delete parameters from the filter state. 
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Figure 1 Filter Subroutine 
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Smoothing 


Optimal smoothing is performed using the backward recursion 
developed by Rauch, Tung and Striebel [4]. The final estimate of the 
filter is used to initialize the smoother equations. The smoother gain 
matrix at time t^ is computed as: 

G k = $ kV r " Q k+l P k+l/k^ 


Then the smoothed state vector and covariance are computed as: 


-k/m ' ?k/k + G k^-k+l/m~ -k+1 /k) 
P k/m = P k/k + G k^ P k+l/m _P k+l/k^ G k 


where the notation 


-i/j 


upon measurements up to time t. 

j 

estimate at time t^ 


means the estimate x at time 
In other words, x,. 


t. based 
is the 


a prtorz 
at time t 


and 


k A k/m 

the last data point). 


k+l/k 

is the a posteriori estimate 
is the smoothed estimate at time t. (t is 


-k/k 


Notice that the gain matrix has the following structure: 

go) G(?r 

_ 0 I _ 

where the partitioning indicated separates the dynamic parameters from 
the biases. Since the number of biases may be several times greater 
than the number of dynamic parameters, the multiplications by 0 or 1 
are avoided in the coding. 
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Although Kalman filter formulations based upon covariance matrices 
are more prone to numerical problems than the factored filters, numerical 
problems are not so severe in the smoother. The smoother equations are 
only evaluated once per mini -batch rather than for each measurement. 
Furthermore, the equations for the smoothed x and P are uncoupled 
since the gain matrix only depends upon variables from the filter. Thus, 
errors in the smoothed P have no effect upon x . 

Disposable Pass Parameters in Smoothing 

It is fairly well known that measurement bias parameters need only 
be included in the filter state during periods when data of the appro- 
priate type is actually being processed. Outside the data interval, the 
solution for the pass parameters has no effect upon the solution for 
the common parameters. 

We are not aware of any published reference which demonstrates that 
the "disposable parameter" approach is also valid for smoothing. 
Therefore, this section shows that the approach is valid and demonstrates 
how it is implemented for the present problem. The following derivation 
is basically the same as that given by Tanenbaum. 1 

Fraser and Potter [2] showed that the optimum smoother could also 
be derived as the linear combination of a forward filter which includes 
a priori information and a backward filter which does not include 
a priori. The results obtained from such a filter will be identical to 
those obtained by the RTS algorithm. 

Consider the case shown in the figure where the forward filter has 
processed data from pass a but not b while the backward filter has 
processed data from b but not a . 


“'Tanenbaum, M. , private communication, NSWC/Dahl gren , December 1977. 
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The filter states and covariances at time t are: 
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where subscript c denotes common parameters. Notice that the a -priori 
information for the pass b parameters of the foward filter is treated 
as if it is a measurement which does not actually enter the forward 
filter until the pass is begun. It can also be shown (with some dif- 
ficulty) that similar results are obtained by allowing it to enter the 
filter at the initial time. The smoothed covariance is obtained as a 
minimum variance combination of the two estimates. Since the errors in 
the two estimates are uncorrelated, the smoothed covariance is simply 
the inverse of the sum of the two information matrices* 


Operations on matrices containing °° must be done with great care. 
The result can, however, be derived more rigorously. 
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P k/m ' < p k 1+p r’)’' 

P cc' P cc +P cc)'’ P cb 
p ac< p cc +p cc>'' P cb 
P ib- P bc< P cc +P cc)'' P cb 

Notice that the solution for the common parameters does not depend 
upon the pass parameters. Furthermore, the solution for pass a does not 
depend upon the pass b parameters (and vice versa). This verifies that 
it is not necessary to carry the pass parameters outside of the pass. 
However, we must also verify that the pass parameters can be "reconstructed" 
in the RTS formulation of the smoother. 


-1 -1 ”1 I 1 

(P >P' ) P' (P +P' ) 1 P 


cc cc 


cc v cc cc' 


ca 


p _p (p +p> ) -1 p 

aa ac v cc cc' ca 
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k/k 


P k+l/k P k+l/k+l 




v»i. 


L k-1 


'k+1 


'k+2 


Consider the case shown in the figure. Assume that the smoothed values 
for t k are to be computed. and P k+1/ , k from the forward fil- 

ter have the same dimension but P k+ iy k+ j does not include the pass 
parameters. Obviously, the smooth covariance, / m » is the same 
dimension as P k+1/k+1 • In the RTS equations, the difference 
P k+l/m " P k+1 'k must corn P u ‘t e d but these two arrays are of different 
dimensions. Therefore, we examine whether the missing terms of aP 
can be reconstructed. Using the results from the forward-backward 
smoother, we find that: 
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AP 


cc 


P k+l/m ' P k+l/k 


ac cc' 


_1 p ) 
cc ca' 

- P cc< P X> 

(p _1 p ) ' 
cc v cc ca' , 

-< p ac p cc> ip cc( p cc p cb’ 


a p = _p (p +p' )" 1 p 

cc cc v cc cc' cc 


where 

partition of P 


is simply computed as the upper left 


k+l/m " r k+l /k 


When written in this form, it is obvious that AP k+ ^ is singular. 
This also shows that the "missing" terms of AP can be reconstructed 
by pre- or post-multiplying by the factor P ac obtained from 
P k+l/k * " Phe rationa ^ for discarding pass parameters after writing the 
filter a priori to the disk should now be obvious. 

By a similar procedure, we can also demonstrate that the pass para- 
meter portion of x k+1 . - x k+1/ , k can be reconstructed as 


A -k+l 


AX_ 

-CC 


( p ac P ci> ‘5cc 


'^k+1 


The equation for the gain matrix requires that P k+ y k be inverted.' 
It can be easily shown [3] that the same results for the smoothed x 
and P will be obtained whether or not the pass parameters are 
included in the gain computation. Thus, the final RTS equations used 
when reconstructing pass parameters are: 


-k/m 

P k/m 


= -k/k 
= p k/k 


+ G ' A -cc 


+ 6' AP 


CC 


g- t 
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where 


G ' = 


r 1 (i-o p" 1 ) 

H CC ^ 'CC CC' 


p P _1 
ac cc 




Examples 

Two examples using simulated data are given to demonstrate the 
improved performance of PREFER. The first is relatively trivial in that 
no modeling errors were included. The test was made simply to evaluate 
the program response to an initial condition error. Table 1 summarizes 
the test case and Figure 2 displays the results. The filter position 
error was initially 20 meters. During the first data pass, the error was 
reduced to 7 meters but during the subsequent data gap, the error rose 
to 38 meters. After the first orbit, the filter error remained below 
1 meter. However, the smoother position error was less than 1.2 meters 
for the enti re run. The smoother error is largest at epoch because 
the 1 sigma a priori error is weighted into the solution. 

The second example is a more rigorous test of the program. It 
includes some additional data types and also has significant force 
modeling errors. Table 2 summarizes the input and Figure 3 displays 
the results. 

The filter estimate has peak errors of 63 meters (mostly cross- 
track) while the maximum error in the smoother estimate is 11.2 meters 
(mostly radial) at the epoch. The peak error in the filter estimate 
occurs at 30 to 40 minutes which corresponds to a minimum error in the 
nominal trajectory. Apparently the filter had an erroneous estimate 
of the gravitational accelerations at the time that a data gap occurred. 
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ORBIT - 350-420 KM ALTITUDE, e = .005, 96.9° INCLINA- 

TION, 180 MINUTES (2 REVOLUTIONS) 

MODEL ERRORS - NONE (NOMINAL TRAJECTORY IS PERFECT) 

TRACKING DATA - 7 GROUND STATIONS, RANGE DATA ONLY, M MEASURE- 

MENT NOISE BUT DATA IS GIVEN A WEIGHT OF 1 METER 

ADJUSTED PARAMETERS - ORBITAL ELEMENTS, MEASUREMENT BIAS AND REFRAC- 
TION PARAMETERS, STATION POSITION ERRORS 

INITIAL CONDITIONS - FILTER ESTIMATE OF SEMI-MAJOR AXIS AT EPOCH IS 

PERTURBED BY 20 METERS do) 


A PRIORI STANDARD DEVIATIONS 

SEMI-MAJOR AXIS - 20 M 

e SIN u - .00001 RADIAN 

e COS u - .00001 RADIAN 

INCLINATION - ,00001 RADIAN 

1 + w - .00001 RADIAN 

n - .00001 RADIAN 

STATION BIAS - 1 M 

STATION REFRACTION - 0,5 M 

STATION POSITION - 5 M (EACH COMPONENT) 

STATE NOISE SPECTRAL DENSITY 

X, Y, Z - .03 M/SEC 1/2 

*, i - .3xl0 _i< M/SEC 372 

Table 1 Summary of Test Case Number 1 
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RSS Position Error (m) 



Figure 2 Error in Estimated Position for Example 1 


ORBIT 


MODEL ERRORS 


TRACKING DATA 


Adjusted 

Parameters 


- 165-264 km altitude, e = ,0075, 96.4° inclination, 
192 MINUTES (2 revolutions) 


- measurement data generated using a 25,25 gravity 

FIELD. NOMINAL TRAJECTORY WAS OBTAINED BY LEAST 
SQUARES FITTING THE TRUE TRAJECTORY USING A 8,8 
GRAVITY FIELD. THE RESULTING POSITION ERRORS ARE 
LESS THAN 53 METERS. ALSO, SINUSOIDAL ERRORS WERE 
ADDED TO THE POSITIONS ON THE GPS TRAJECTORY FILE. 
THE STANDARD DEVIATIONS FOR THE PEAK ERRORS WERE: 
10 METERS ALONG-TRACK, 6 METERS CROSS-TRACK AND 
2 METERS RADIALLY. 


- 16 GROUND STATIONS: ALL HAVE RANGE DATA BUT TWO 

•ALSO HAVE RANGE DIFFERENCE AND ANOTHER TWO HAVE 
DOPPLER DATA. DATA IS NOISELESS BUT IS GIVEN 
WEIGHTS OF 1 METER (RANGE), 6 CM (RANGE DIFFERENCE) 
AND 0.2x10"^ (DOPPLER). 6 GPS SATELLITES (pseudo 
RANGE AND DELTA-RANGE). DATA HAS MEASUREMENT NOISE 
OF 1.5 METERS (PSEUDO RANGE) AND 2 CM (PSEUDO DELTA- 
RANGE). DATA IS WEIGHTED ACCORDINGLY, 


^Gravitational acceleration, host clock errors, 

STATION MEASUREMENT BIASES AND REFRACTION, STATION * 
POSITIONS, GPS POSITIONS AND TIMING, 


Table 2 Summary of Test Case Number 2 



RSS Position Error (m) 



Figure 3 Error in Estimated Position for Example 2 



Thus, the error quickly increased until more tracking was obtained. 
Hov/ever, the filter covariance matrix during the data gaps was also 
large so that the smoother could correctly weight the filter estimates. 

Notice that both the filter and smoother estimates are quite 
accurate during the periods when GPS tracking is available. During 
these periods, the smoother estimation error was generally less than 
three meters and the radial component was accurate to within 1.5 meters 
Even during the data gaps, the smoother radial error did not exceed 6 
meters (except at the epoch). This large error occurred at 102 minutes 
from epoch and the nominal trajectory at this time had a 50 meter cross 
track error. 

It should be noted that no great attempt was made to "fine tune" 
the input parameters for this example. Presumably the errors could be 
reduced further by the appropriate choice of state noise variances, 
time constants, etc. 

Summary 

The results of the various tests on simulated data demonstrate 
that PREFER has great potential for improving orbit determination of 
low altitude satellites. 
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